find_new_clusters

Parameters

suffix = ""
data_to_read = "./Data/acc_tpm_nCount_mito_no146_15k_cancercells.rds"

functions

Data

acc = readRDS(file = data_to_read)
(message("reading '" %>% paste0(data_to_read %>% basename()) %>% paste0("'")))

reading 'acc_tpm_nCount_mito_no146_15k_cancercells.rds'

NULL

pathways_scores = fread(file = "./Data/ACC_Canonical_Pathway_Scores.txt",sep = ",") %>% as.matrix(rownames=1) %>% t() %>%  as.data.frame()
hallmark_scores = fread(file = "./Data/ACC_Hallmark_Pathway_Scores.txt",sep = ",") %>% as.matrix(rownames=1) %>% t() %>%  as.data.frame()

ln_list = c("ACC22.LN.P11", "ACC22.P12.LN","ACC7.P13")
ln_plates = FetchData(object = acc,vars = "orig.ident") %>% mutate(
  tumor_type = if_else(condition = orig.ident %in% ln_list
                       ,true = "LN"
                       ,false = "primary"))

ln_plates["orig.ident"] <-NULL
acc= AddMetaData(object = acc,metadata = ln_plates)

pathways_scores = cbind(pathways_scores,hallmark_scores)
pathways_scores = pathways_scores[ , colSums(is.na(pathways_scores))==0] #remove cols with NA
pathways_scores = pathways_scores [rownames(pathways_scores) %in% colnames(acc),] #remove cells not in dataset
pathways_scores =  pathways_scores[order(row.names(pathways_scores)),] #order cells like dataset

# run-dim-reduction on genes:
acc <- FindVariableFeatures(acc, selection.method = "vst", nfeatures = 2000)

acc <- ScaleData(acc)

acc <- RunPCA(acc)

ElbowPlot(acc)

acc <- RunUMAP(acc, dims = 1:5)

pathway_scores_assay <- CreateAssayObject(counts = pathways_scores %>% t()) #create an assay

Warning: Feature names cannot have underscores ('_'), replacing with dashes ('-')

acc[["pathway_scores"]] = pathway_scores_assay

Warning: Keys should be one or more alphanumeric characters followed by an underscore, setting key from pathway_scores_ to pathwayscores_

acc <- RunUMAP(acc, dims = 1:5,reduction ="PCA_pathway_scores",reduction.name = "pathway_scores_umap")

acc umaps

gene expression

DimPlot(acc,group.by = "patient.ident")

pathways scores

DimPlot(acc,reduction = "pathway_scores_umap",group.by = "patient.ident")

UMAP by tumor type

genes

DimPlot(acc,reduction = "umap",group.by = "tumor_type")

pathwyas

DimPlot(acc,reduction = "pathway_scores_umap",group.by = "tumor_type")

All PC’s

for (i in 1:8) {
  cat("### PC",i," \n")
  print(VizDimLoadings(acc, dims = i, reduction = "PCA_pathway_scores"))
  cat(' \n\n')
  plot.new()
  dev.off()
}

PC 1

PC 2

PC 3

PC 4

PC 5

PC 6

PC 7

PC 8

cycling cells clustring

hallmark_name = "HALLMARK_G2M_CHECKPOINT"
genesets  =GSEABase::getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
var_features=acc@assays$RNA@var.features
geneIds= genesets[[hallmark_name]]@geneIds
score <- apply(acc@assays$RNA@data[intersect(geneIds,var_features),],2,mean)
acc=AddMetaData(acc,score,hallmark_name)

genes

FeaturePlot(acc, reduction = "umap",features = "HALLMARK_G2M_CHECKPOINT")

pathwyas

FeaturePlot(acc, reduction = "pathway_scores_umap",features = "HALLMARK_G2M_CHECKPOINT")

UMAP clusters

acc <- FindNeighbors(acc, dims = 1:10,reduction = "PCA_pathway_scores")
acc <- FindClusters(acc, resolution = 0.1,graph.name = "pathway_scores_snn")

DimPlot(acc,reduction = "pathway_scores_umap")

DEG

markers = FindMarkers(object = acc,ident.1 = "0",ident.2 = "1",assay = "pathway_scores",min.cells.feature = 10,logfc.threshold = 0)

markers

Top 10 pathways

DefaultAssay(acc )= "pathway_scores"
for (i in 1:10) {
  pathway = rownames(markers)[i]
  cat("### ",i," \n")
  print(
    FeaturePlot(object = acc,features = pathway,reduction = "pathway_scores_umap")
    )
  plot.new()
  dev.off()
  cat(' \n\n')
}

1

2

3

4

5

6

7

8

9

10

NA

HEAD-AND-NECK-SQUAMOUS

FeaturePlot(object = acc,features = "WP-HEAD-AND-NECK-SQUAMOUS-CELL-CARCINOMA",reduction = "pathway_scores_umap")

gs=acc@assays$RNA@var.features

myoscore=apply(acc@assays$RNA@scale.data[intersect(c("TP63","TP73","CAV1","CDH3","KRT5","KRT14","ACTA2","TAGLN","MYLK","DKK3"),gs),],2,mean)

lescore=apply(acc@assays$RNA@scale.data[intersect(c("KIT","EHF","ELF5","KRT7","CLDN3","CLDN4","CD24","LGALS3","LCN2","SLPI"),gs),],2,mean)
acc=AddMetaData(acc,lescore-myoscore,"luminal_over_myo")

#set lum_or_myo metadata
luminal_over_myo = FetchData(object = acc,vars = "luminal_over_myo")
luminal_over_myo$lum_or_myo  = case_when(luminal_over_myo$luminal_over_myo >1~"lum",luminal_over_myo$luminal_over_myo <(-1)~"myo",TRUE~"NA")
luminal_over_myo$luminal_over_myo <-NULL
acc=AddMetaData(object = acc,metadata = luminal_over_myo,col.name = "lum_or_myo")

FeaturePlot(object = acc,features = "luminal_over_myo",reduction = "pathway_scores_umap")

DimPlot(acc,group.by = "lum_or_myo",cols = c("red","green","grey"))